Accreditation exercise for 16S rRNA data analysis

Pre-processing

Cleaning the taxonomy and feature abundance file.

sed 's/#OTU//' feature-table.tsv > table.tsv

asvs_abundance <- read.table("table.tsv", header = TRUE, row.names = "ID")
# dim(asvs_abundance)

sed 's/\[//g' taxonomy.tsv | sed 's/\]//g' > tax.tsv

asv_taxonomy <- read.delim("tax.tsv", row.names = 1)

asv_taxonomy[] <- lapply(asv_taxonomy, function(x) gsub("D_0__|D_1__|D_2__|D_3__|D_4__|D_5__|D_6__|D_7__|D_8__|D_9__|D_10__|D_11__|D_12__|D_13__|D_14__", 
    "", x))
asv_taxonomy <- asv_taxonomy[asv_taxonomy$Confidence > 0.97, ]  #consider 97% classification confidence threshold
# Re-organise taxonomy dataframe to have seven columns each with single
# taxonomic rank
require(splitstackshape)
tmp <- cSplit(asv_taxonomy, "Taxon", ";")
tmp <- tmp[, -1]  #remove confidence column
tmp <- tmp[, 1:7]  #pick first seven columns
names(tmp) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")  #assign new column names
tax <- data.frame(tmp, row.names = rownames(asv_taxonomy))
tax <- tax[tax$Phylum != "" & !is.na(tax$Phylum), ]  #remove ASVs without Phylum level classification
tax$Species[is.na(tax$Species)] <- ""
tax$Genus[is.na(tax$Genus)] <- ""
tax$G_species <- paste(tax$Genus, tax$Species, sep = " ")  #create a new taxonomic rank (species and genus for which species is not assigned)

Load phylogenetic tree

require(ape)
asv_tree <- read_tree("tree.nwk")

Construct a phyloseq object, composed of ASV abundance, taxonomy, tree and meta data.

require(phyloseq)
ASV <- otu_table(as.matrix(asvs_abundance), taxa_are_rows = T)
TAX <- tax_table(as.matrix(tax))
meta_data <- read.delim("set2_meta.txt", row.names = 1)
meta_data$sample <- row.names(meta_data)
meta_data$BV <- as.factor(meta_data$BV)
meta_data$Inflammation <- as.factor(meta_data$Inflammation)
SAM <- sample_data(meta_data)
tree <- compute.brlen(asv_tree, method = "Grafen")
physeq <- merge_phyloseq(phyloseq(ASV, TAX), SAM, tree)

Obtain number and percentage of ASVs identified at each taxonomic rank

sp <- dim(tax[!is.na(tax$Species), ])[1]  #Species level
sp/dim(tax)[1]
gp <- dim(tax[!is.na(tax$Genus) & tax$Genus != "", ])[1]  #Genus level
gp/dim(tax)[1]
fm <- dim(tax[!is.na(tax$Family) & tax$Family != "", ])[1]  #Family level
fm/dim(tax)[1]
saveRDS(physeq, "physeq_silva_972")  #save the phyloseq object for future use

From now onwards, we use the previous saved phyloseq object for downstream analysis. The phyloseq object used downstream is available at https://github.com/galaxyuvri-ea/16S-rRNA-analysis/raw/master/dataFiles/physeq_silva_972.

require(kableExtra)
require(phyloseq)
phy <- readRDS("physeq_silva_972")
meta <- data.frame(sample_data(phy))
meta %>% kable(booktabs = TRUE, longtable = TRUE) %>% kable_styling(latex_options = c("hold_position", 
    "repeat_header", position = "")) %>% scroll_box(width = "100%", height = "300px")

require(kableExtra)
stats <- read.table("data/stats.tsv", header = T)
stats %>% kable(booktabs = TRUE, longtable = TRUE) %>% kable_styling(latex_options = c("hold_position", 
    "repeat_header", position = "")) %>% scroll_box(width = "100%", height = "300px")

mean(stats$input)

[1] 185871.8

mean(stats$filtered)

[1] 143041.2

mean(stats$merged)

[1] 137931.7

mean(stats$non.chimeric)

[1] 137931.7

lost_filt <- (mean(stats$input) - mean(stats$filtered))/mean(stats$input)
lost_filt

[1] 0.2304308

lost_merge <- (mean(stats$filtered) - mean(stats$merged))/mean(stats$filtered)
lost_merge

[1] 0.03572084

lost_chime <- (mean(stats$merged) - mean(stats$non.chimeric))/mean(stats$merged)
lost_chime

[1] 0

lost_overall <- (mean(stats$input) - mean(stats$non.chimeric))/mean(stats$input)
lost_overall

[1] 0.2579204

require(plotly)
df <- stats[, -c(1)]  #sample ids
p <- ggplot(stack(df), aes(x = ind, y = values)) + geom_boxplot(col = "blue", 
    outline = FALSE) + theme_bw() + ylab("Number of reads") + xlab("")
ggplotly(p)

df <- as.data.frame(sample_sums(phy))
colnames(df) <- c("Number_of_ASVs")
df %>% kable(booktabs = TRUE, longtable = TRUE) %>% kable_styling(latex_options = c("hold_position", 
    "repeat_header", position = "")) %>% scroll_box(width = "100%", height = "200px")

phy <- subset_samples(phy, !is.na(Inflammation))  #remove samples for which without values for inflammation
phy <- subset_samples(phy, sample_sums(phy) > 5000)  #retain samples with 5000 and more reads.
phy <- subset_samples(phy, sample != "S120")  #remove outlier sample 
phy

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3222 taxa and 111 samples ]
sample_data() Sample Data:       [ 111 samples by 5 sample variables ]
tax_table()   Taxonomy Table:    [ 3222 taxa by 8 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 3222 tips and 3218 internal nodes ]

total <- median(sample_sums(phy))
standf <- function(x, t = total) round(t * (x/sum(x)))
M.std <- transform_sample_counts(phy, standf)

M.std <- filter_taxa(M.std, function(x) sum(x > 10) > (0.01 * length(x)) | sum(x) > 
    0.001 * total, TRUE)
ntaxa(M.std)

[1] 1075

Unsupervised clustering

require(cluster)
library(factoextra)
df1 <- data.frame(otu_table(M.std))
df_viz <- fviz_nbclust(df1, kmeans, method = "wss") + theme_minimal()
res_fanny <- fanny(t(df1), k = 3, metric = "SqEuclidean")
sample_data(M.std)$cluster <- as.character(res_fanny$clustering)

require(microbiomeSeq)
require(pheatmap)
phyTop <- taxa_level(M.std, "G_species")  #reduce dataset to species and Genus rank
TopASVs <- names(sort(taxa_sums(phyTop), TRUE))[1:50]  #pick 50 most abundant taxa
df1 <- as.data.frame(otu_table(t(phyTop)))
df1 <- df1[which(rownames(df1) %in% TopASVs), ]
df1 <- df1[which(rownames(df1) != " "), ]
df <- data.frame((sample_data(phyTop))[, c("Inflammation", "BV")])
p <- pheatmap(log2(df1 + 1), cluster_rows = FALSE, show_rownames = TRUE, show_colnames = F, 
    cluster_cols = TRUE, annotation_col = df, annotation_row = NA)
p

Alpha diversity

require(plotly)
p <- plot_richness(M.std, x = "cluster", color = "cluster", measures = c("Observed", 
    "Simpson"), title = paste0("Standardized to total reads"))
p <- p + geom_boxplot() + geom_point() + theme_bw() + theme(axis.text = element_text(size = 10, 
    face = "bold"), axis.title = element_text(size = 16, face = "bold")) + geom_point(size = 2)
ggplotly(p)

p <- plot_richness(M.std, x = "Inflammation", color = "Inflammation", measures = c("Observed", 
    "Simpson"), title = paste0("Standardized to total reads"))
p <- p + geom_boxplot() + geom_point() + theme_bw() + theme(axis.text = element_text(size = 10, 
    face = "bold"), axis.title = element_text(size = 16, face = "bold")) + geom_point(size = 2)
ggplotly(p)

p <- plot_richness(M.std, x = "BV", color = "BV", measures = c("Observed", "Simpson"), 
    title = paste0("Standardized to total reads"))
p <- p + geom_boxplot() + geom_point() + theme_bw() + theme(axis.text = element_text(size = 10, 
    face = "bold"), axis.title = element_text(size = 16, face = "bold")) + geom_point(size = 2)
ggplotly(p)

est <- estimate_richness(M.std, split = TRUE, measures = c("Simpson"))
BV_div <- cbind(est, sample_data(M.std)[, "BV"])
t <- wilcox.test(BV_div$Simpson ~ BV_div$BV)
t

    Wilcoxon rank sum test with continuity correction

data:  BV_div$Simpson by BV_div$BV
W = 398, p-value = 1.71e-11
alternative hypothesis: true location shift is not equal to 0

BInflammation_div <- cbind(est, sample_data(M.std)[, "Inflammation"])
t <- wilcox.test(BInflammation_div$Simpson ~ BInflammation_div$Inflammation)
t

    Wilcoxon rank sum test with continuity correction

data:  BInflammation_div$Simpson by BInflammation_div$Inflammation
W = 2204, p-value = 3.854e-05
alternative hypothesis: true location shift is not equal to 0

Cluster_div <- cbind(est, sample_data(M.std)[, "cluster"])
dunn.test::dunn.test(Cluster_div$Simpson, Cluster_div$cluster)

  Kruskal-Wallis rank sum test

data: x and group
Kruskal-Wallis chi-squared = 58.3445, df = 2, p-value = 0

                           Comparison of x by group                            
                                (No adjustment)                                
Col Mean-|
Row Mean |          1          2
---------+----------------------
       2 |  -5.906293
         |    0.0000*
         |
       3 |  -0.080297   6.501771
         |     0.4680    0.0000*

alpha = 0.05
Reject Ho if p <= alpha/2

Beta diversity

# ordination based on NMDS and bray-curtis distance
NMDS_bray <- ordinate(M.std, "NMDS")

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2055082 
Run 1 stress 0.2309574 
Run 2 stress 0.217807 
Run 3 stress 0.2221731 
Run 4 stress 0.1968815 
... New best solution
... Procrustes: rmse 0.0485541  max resid 0.2755098 
Run 5 stress 0.222773 
Run 6 stress 0.2169525 
Run 7 stress 0.2262007 
Run 8 stress 0.226073 
Run 9 stress 0.2171289 
Run 10 stress 0.2318284 
Run 11 stress 0.2278323 
Run 12 stress 0.2134497 
Run 13 stress 0.2192831 
Run 14 stress 0.2118285 
Run 15 stress 0.2277937 
Run 16 stress 0.2111281 
Run 17 stress 0.2105766 
Run 18 stress 0.4130426 
Run 19 stress 0.2136909 
Run 20 stress 0.2176685 
*** No convergence -- monoMDS stopping criteria:
    20: stress ratio > sratmax

title = c("NMDS of 16S microbiome, Bray-Curtis distance")

NMDS.bray.plot <- plot_ordination(M.std, NMDS_bray, color = "BV", shape = "Inflammation", 
    title = title)
NMDS.bray.plot <- NMDS.bray.plot + theme(axis.text = element_text(size = 16, 
    face = "bold"), axis.title = element_text(size = 18, face = "bold"), legend.title = element_text(size = 14)) + 
    theme_bw() + labs(color = "BV") + geom_point(size = 5)
NMDS.bray.plot

PCoA.ord.bray <- ordinate(M.std, "PCoA", "bray")
title = c("PCoA of 16S microbiome, Bray-Curtis distance")
PCoA.ord <- plot_ordination(M.std, PCoA.ord.bray, color = "cluster", shape = "BV", 
    title = title)
PCoA.ord <- PCoA.ord + theme(axis.text = element_text(size = 16, face = "bold"), 
    axis.title = element_text(size = 18, face = "bold"), legend.title = element_text(size = 14)) + 
    theme_bw() + labs(color = "cluster") + geom_point(size = 5)
PCoA.ord

# permanova
require(vegan)
phy_bray <- phyloseq::distance(M.std, method = "bray")
sampledf <- data.frame(sample_data(M.std))

adonis(phy_bray ~ BV, data = sampledf)

Call:
adonis(formula = phy_bray ~ BV, data = sampledf) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
BV          1     7.396  7.3956  27.198 0.19969  0.001 ***
Residuals 109    29.639  0.2719         0.80031           
Total     110    37.035                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

beta <- betadisper(phy_bray, sampledf$BV)
permutest(beta)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups      1 0.00977 0.0097732 0.5461    999  0.469
Residuals 109 1.95083 0.0178975                     

adonis(phy_bray ~ Inflammation, data = sampledf)

Call:
adonis(formula = phy_bray ~ Inflammation, data = sampledf) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
Inflammation   1     3.676  3.6759  12.011 0.09926  0.001 ***
Residuals    109    33.359  0.3060         0.90074           
Total        110    37.035                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

beta <- betadisper(phy_bray, sampledf$Inflammation)
permutest(beta)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups      1 0.00025 0.0002487 0.0177    999  0.904
Residuals 109 1.53154 0.0140509                     

adonis(phy_bray ~ cluster, data = sampledf)

Call:
adonis(formula = phy_bray ~ cluster, data = sampledf) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
cluster     2    17.856  8.9278  50.273 0.48213  0.001 ***
Residuals 108    19.179  0.1776         0.51787           
Total     110    37.035                 1.00000           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

beta <- betadisper(phy_bray, sampledf$BV)
permutest(beta)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups      1 0.00977 0.0097732 0.5461    999  0.448
Residuals 109 1.95083 0.0178975                     

Multivariate logistic regression

sampledf$cluster <- as.factor(sampledf$cluster)
logit_model <- glm(Inflammation ~ cluster + BMI + Age + BV, family = binomial(link = "logit"), 
    data = sampledf)
summary(logit_model)

Call:
glm(formula = Inflammation ~ cluster + BMI + Age + BV, family = binomial(link = "logit"), 
    data = sampledf)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8491  -0.6003  -0.5532   0.8470   1.9964  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  1.33330    3.24810   0.410  0.68145   
cluster2    -0.85151    0.88297  -0.964  0.33486   
cluster3    -0.61949    0.68378  -0.906  0.36495   
BMI          0.02086    0.05698   0.366  0.71434   
Age         -0.02631    0.15561  -0.169  0.86574   
BVPositive  -2.22065    0.76632  -2.898  0.00376 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.19  on 99  degrees of freedom
Residual deviance: 101.90  on 94  degrees of freedom
  (11 observations deleted due to missingness)
AIC: 113.9

Number of Fisher Scoring iterations: 4

exp(logit_model$coefficients)

(Intercept)    cluster2    cluster3         BMI         Age  BVPositive 
  3.7935318   0.4267703   0.5382209   1.0210768   0.9740343   0.1085389 

exp(confint(logit_model))

                  2.5 %       97.5 %
(Intercept) 0.005649307 2196.3498153
cluster2    0.075719698    2.6415210
cluster3    0.132252489    2.0069491
BMI         0.912317941    1.1434501
Age         0.717828619    1.3306259
BVPositive  0.020596434    0.4536472

Wald test:
----------

Chi-squared test:
X2 = 1.1, df = 2, P(> X2) = 0.57

Differential abundance testing

# Differential abundance using DESEQ2
require(DESeq2)
M.f.des <- taxa_level(M.std, "G_species")
dds <- phyloseq_to_deseq2(M.f.des, ~Inflammation)
# geometric mean, set to zero when all coordinates are zero
geo_mean_protected <- function(x) {
    if (all(x == 0)) {
        return(0)
    }
    exp(mean(log(x[x != 0])))
}
geoMeans <- apply(counts(dds), 1, geo_mean_protected)
dds <- estimateSizeFactors(dds, geoMeans = geoMeans)
dds <- DESeq(dds, fitType = "local")

res = results(dds, cooksCutoff = FALSE)
df <- as.data.frame(colData(dds)[, "BV"])
colnames(df) <- "BV"

alpha = 0.01
sigtab = as.data.frame(res[which(res$padj < alpha & abs(res$log2FoldChange) > 
    2), ])
# sigtab = cbind(as(sigtab, 'data.frame'),
# as(tax_table(M.f.des)[rownames(sigtab), ], 'matrix'))
dim(sigtab)

[1] 14  6

table_sig <- sigtab[c("log2FoldChange", "padj")]
table_sig$Abundant_Group <- levels(df$BV)[as.numeric(sigtab$log2FoldChange > 
    0) + 1]
table_sig %>% kable(booktabs = TRUE, longtable = TRUE) %>% kable_styling(latex_options = c("hold_position", 
    "repeat_header", position = "")) %>% scroll_box(width = "100%", height = "300px")

library("pheatmap")
select <- order(rowMeans(counts(dds, normalized = TRUE)), decreasing = TRUE)[1:50]
df <- as.data.frame(colData(dds)[, c("Inflammation", "BV", "cluster")])

df1 <- data.frame(assay(dds))[select, ]
p <- pheatmap(log2(df1 + 1), cluster_rows = FALSE, show_rownames = TRUE, show_colnames = F, 
    cluster_cols = TRUE, annotation_col = df, annotation_row = NA)
p

require(randomForest)
data.rf <- data.frame(otu_table(M.f.des))
data.rf$BV <- sample_data(M.std)$BV
BV.rf <- randomForest(BV ~ ., data = data.rf, importance = TRUE, proximity = TRUE)
df_accuracy <- as.data.frame(BV.rf$importance)
df_accuracy$Taxa <- rownames(df_accuracy)
df_accuracy <- df_accuracy[order(df_accuracy$MeanDecreaseAccuracy, decreasing = TRUE), 
    ][1:20, ]
df_accuracy$Taxa <- factor(df_accuracy$Taxa, levels = df_accuracy$Taxa)
mda_plot <- ggplot(data = df_accuracy, aes(x = Taxa, y = MeanDecreaseAccuracy)) + 
    theme_bw()
mda_plot <- mda_plot + geom_bar(stat = "identity", fill = "darkblue", width = 0.5)
mda_plot <- mda_plot + theme(axis.text.x = element_text(angle = 90, hjust = 1, 
    vjust = 0.5))
mda_plot <- mda_plot + xlab("") + ylab("Mean Decrease in Accuracy") + coord_flip() + 
    theme(axis.text = element_text(size = 10, face = "bold"), axis.title = element_text(size = 16, 
        face = "bold"))
mda_plot

Session Info

R version 3.4.1 (2017-06-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] randomForest_4.6-14        DESeq2_1.16.1             
 [3] SummarizedExperiment_1.6.5 DelayedArray_0.2.7        
 [5] matrixStats_0.54.0         Biobase_2.36.2            
 [7] GenomicRanges_1.28.6       GenomeInfoDb_1.12.3       
 [9] IRanges_2.10.5             S4Vectors_0.14.7          
[11] BiocGenerics_0.22.1        aod_1.3.1                 
[13] vegan_2.5-4                lattice_0.20-38           
[15] permute_0.9-5              pheatmap_1.0.12           
[17] microbiomeSeq_0.1          factoextra_1.0.5          
[19] cluster_2.0.7-1            plotly_4.9.0              
[21] ggplot2_3.2.0              phyloseq_1.25.2           
[23] kableExtra_1.1.0           rmdformats_0.3.5          
[25] knitr_1.23                

loaded via a namespace (and not attached):
 [1] backports_1.1.4      uuid_0.1-2           Hmisc_4.2-0         
 [4] plyr_1.8.4           igraph_1.2.4         lazyeval_0.2.2      
 [7] sp_1.3-1             splines_3.4.1        crosstalk_1.0.0     
[10] BiocParallel_1.10.1  rncl_0.8.3           digest_0.6.19       
[13] foreach_1.4.4        htmltools_0.3.6      gdata_2.18.0        
[16] memoise_1.1.0        checkmate_1.9.3      magrittr_1.5        
[19] annotate_1.54.0      Biostrings_2.44.2    readr_1.3.1         
[22] gmodels_2.18.1       prettyunits_1.0.2    colorspace_1.4-1    
[25] blob_1.1.1           rvest_0.3.4          ggrepel_0.8.1       
[28] xfun_0.8             dplyr_0.8.0.1        crayon_1.3.4        
[31] RCurl_1.95-4.12      jsonlite_1.6         genefilter_1.58.1   
[34] phylobase_0.8.6      survival_2.44-1.1    iterators_1.0.10    
[37] ape_5.3              glue_1.3.1           gtable_0.3.0        
[40] zlibbioc_1.22.0      XVector_0.16.0       webshot_0.5.1       
[43] seqinr_3.4-5         questionr_0.7.0      dunn.test_1.3.5     
[46] adegraphics_1.0-15   scales_1.0.0         DBI_1.0.0           
[49] miniUI_0.1.1.1       Rcpp_1.0.1           htmlTable_1.13.1    
[52] viridisLite_0.3.0    xtable_1.8-4         progress_1.2.2      
[55] spData_0.3.0         bit_1.1-14           foreign_0.8-71      
[58] spdep_0.7-9          Formula_1.2-3        htmlwidgets_1.3     
[61] httr_1.4.0           RColorBrewer_1.1-2   acepack_1.4.1       
[64] pkgconfig_2.0.2      XML_3.98-1.20        nnet_7.3-12         
[67] deldir_0.1-16        locfit_1.5-9.1       labeling_0.3        
[70] AnnotationDbi_1.38.2 tidyselect_0.2.5     rlang_0.3.1         
[73] reshape2_1.4.3       later_0.8.0          munsell_0.5.0       
 [ reached getOption("max.print") -- omitted 56 entries ]

Acknowledgement

• H3ABioNet
• Uganda Virus Research Institute (UVRI), BCB Hub
• UCT-CBiO

References